欢迎关注“小丫画图”公众号，回复“小白”，看小视频，实现点鼠标跑代码。

小丫微信: epigenomics E-mail:

作者：大鱼海棠，他的更多作品看这里https://k.koudai.com/OFad8N0w

单位：中国药科大学生物统计和计算药学研究中心，国家天然药物重点实验室

小丫编辑校验

需求描述

主要是红框中的图，都是公共数据和已发表的文章

We then compared our classification based on multi-platform data from the TMDU and TCGA studies

根据以往的分子分型marker预测现有样本属于哪一类，再跟当前分型比较，计算p值。

出自https://linkinghub.elsevier.com/retrieve/pii/S2352396418306340

Fig. 3. Summary ofmolecular classification of HCC. (a) Comparison of aggregate scores with gene sets associated with the previously definedmolecular classifications of HCC in the TMDU test study (upper) and TCGA validation study (lower). (b) Schematic representation of molecular subtypes.

2017年的这篇Cell文章里就出现了类似的图：

出自https://linkinghub.elsevier.com/retrieve/pii/S0092867417306396

Figure 2. Liver Cancers Show Distinct Gene Hypermethylation Patterns (A) Unsupervised clustering analysis of gene hypermethylation in HCC tumor relative to normal tissue reveals four distinct subgroups. Roughly 15,000 CpG sites showing significant hypermethylation in 196 HCC patients were analyzed and are shown in heatmap format with normal tissues and tumors organized in columns according to cluster designation. Intensity of methylation for each CpG site is indicated by row. Above the heatmap the four distinct hypermethylation clusters are shown, and below are bars indicating the distribution of clinical and molecular attributes of the individual tumors by cluster. To the right, p values indicate significant non-random distributions for each attribute.

应用场景

根据表达谱对HCC样本进行聚类，并和其他已有的分类器做比较。

环境设置

使用国内镜像安装包

options("repos"= c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
options(BioC_mirror="http://mirrors.tuna.tsinghua.edu.cn/bioconductor/")
devtools::install_github("Lothelab/CMScaller")

加载包

library(CMScaller) # 用来实现nearest template prediction(NTP)
## CMScaller v2.0.1; genome annotation: GENCODE v32/GRCh38.p13
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.1.2
## ── Attaching packages
## ───────────────────────────────────────
## tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6     ✔ purrr   0.3.4
## ✔ tibble  3.1.8     ✔ dplyr   1.0.9
## ✔ tidyr   1.2.0     ✔ stringr 1.4.0
## ✔ readr   2.1.2     ✔ forcats 0.5.1
## Warning: package 'ggplot2' was built under R version 4.1.2
## Warning: package 'tibble' was built under R version 4.1.2
## Warning: package 'tidyr' was built under R version 4.1.2
## Warning: package 'readr' was built under R version 4.1.2
## Warning: package 'dplyr' was built under R version 4.1.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(clusterProfiler)
## 
## Registered S3 method overwritten by 'ggtree':
##   method      from 
##   identify.gg ggfun
## clusterProfiler v4.0.5  For help: https://yulab-smu.top/biomedical-knowledge-mining-book/
## 
## If you use clusterProfiler in published research, please cite:
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation. 2021, 2(3):100141. doi: 10.1016/j.xinn.2021.100141
## 
## Attaching package: 'clusterProfiler'
## 
## The following object is masked from 'package:purrr':
## 
##     simplify
## 
## The following object is masked from 'package:stats':
## 
##     filter
library(org.Hs.eg.db)
## Loading required package: AnnotationDbi
## Loading required package: stats4
## Loading required package: BiocGenerics
## 
## Attaching package: 'BiocGenerics'
## 
## The following objects are masked from 'package:dplyr':
## 
##     combine, intersect, setdiff, union
## 
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## 
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, basename, cbind, colnames,
##     dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
##     grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
##     rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
##     union, unique, unsplit, which.max, which.min
## 
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## 
## Loading required package: IRanges
## Loading required package: S4Vectors
## 
## Attaching package: 'S4Vectors'
## 
## The following object is masked from 'package:clusterProfiler':
## 
##     rename
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, rename
## 
## The following object is masked from 'package:tidyr':
## 
##     expand
## 
## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname
## 
## 
## Attaching package: 'IRanges'
## 
## The following object is masked from 'package:clusterProfiler':
## 
##     slice
## 
## The following objects are masked from 'package:dplyr':
## 
##     collapse, desc, slice
## 
## The following object is masked from 'package:purrr':
## 
##     reduce
## 
## 
## Attaching package: 'AnnotationDbi'
## 
## The following object is masked from 'package:clusterProfiler':
## 
##     select
## 
## The following object is masked from 'package:dplyr':
## 
##     select
library(ClassDiscovery)
## Loading required package: cluster
## Warning: package 'cluster' was built under R version 4.1.2
## Loading required package: oompaBase
library(ComplexHeatmap)
## Loading required package: grid
## ========================================
## ComplexHeatmap version 2.8.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
## 
## If you use it in published research, please cite:
## Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
##   genomic data. Bioinformatics 2016.
## 
## The new InteractiveComplexHeatmap package can directly export static 
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
## 
## This message can be suppressed by:
##   suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
library(gplots)
## Warning: package 'gplots' was built under R version 4.1.2
## 
## Attaching package: 'gplots'
## 
## The following object is masked from 'package:oompaBase':
## 
##     redgreen
## 
## The following object is masked from 'package:IRanges':
## 
##     space
## 
## The following object is masked from 'package:S4Vectors':
## 
##     space
## 
## The following object is masked from 'package:stats':
## 
##     lowess
Sys.setenv(LANGUAGE = "en") #显示英文报错信息
options(stringsAsFactors = FALSE) #禁止chr转成factor

自定义函数

standarize.fun <- function(indata=NULL, halfwidth=NULL, centerFlag=T, scaleFlag=T) {  
  outdata=t(scale(t(indata), center=centerFlag, scale=scaleFlag))
  if (!is.null(halfwidth)) {
    outdata[outdata>halfwidth]=halfwidth
    outdata[outdata<(-halfwidth)]= -halfwidth
  }
  return(outdata)
}

fpkmToTpm <- function(fpkm)
{
  exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}

输入文件

data_mutations.txt，突变数据，下载自cBioPortalhttps://www.cbioportal.org/。

LIHC.htseq_fpkm.tsv.gz，表达谱数据FPKM，已经过log2(fpkm+1)转换，下载地址：https://xenabrowser.net/datapages/?dataset=TCGA-LIHC.htseq_fpkm.tsv&host=https%3A%2F%2Fgdc.xenahubs.net&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443

gencode.v22.annotation.gene.probeMap，ID/Gene Mapping，下载地址同上。

# 读取突变数据
maf <- read_tsv("data_mutations.txt", comment = "#")
## Rows: 223 Columns: 114
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (100): Hugo_Symbol, Center, NCBI_Build, Strand, Consequence, Variant_Cla...
## dbl  (14): Entrez_Gene_Id, Chromosome, Start_Position, End_Position, t_ref_c...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# 把突变数据转成01矩阵，方法跟FigureYa288MutualExclusivity一样
mut.binary <- matrix(0,nrow = length(unique(maf$Hugo_Symbol)),ncol = length(unique(maf$Tumor_Sample_Barcode)),dimnames = list(unique(maf$Hugo_Symbol),unique(maf$Tumor_Sample_Barcode)))
for (i in colnames(mut.binary)) {
  tmp <- maf[which(maf$Tumor_Sample_Barcode == i),]
  tmp <- tmp[which(tmp$Variant_Classification %in% c("Frame_Shift_Del", "Frame_Shift_Ins", "Splice_Site", "Translation_Start_Site","Nonsense_Mutation", "Nonstop_Mutation", "In_Frame_Del","In_Frame_Ins", "Missense_Mutation")),]
  for (j in tmp$Hugo_Symbol)
    mut.binary[j,i] <- 1
}
mut.binary <- as.data.frame(mut.binary)

# 读取表达谱
fpkm <- read.delim("TCGA-LIHC.htseq_fpkm.tsv.gz", sep = "\t", row.names = 1, check.names = F, stringsAsFactors = F, header = T)
fpkm <- 2^fpkm - 1
# 把fpkms转为tpm
tpm <- apply(fpkm, 2, fpkmToTpm)
tpm <- as.data.frame(log2(tpm + 1))
rm(fpkm); gc()
##            used  (Mb) gc trigger   (Mb) limit (Mb)  max used   (Mb)
## Ncells  6028474 322.0   17595237  939.7         NA  15886681  848.5
## Vcells 44534313 339.8  201092013 1534.3      16384 251358087 1917.8
# 加载基因注释并对重复基因名取中位数
Ginfo <- read.delim("gencode.v22.annotation.gene.probeMap",row.names = 1,sep = "\t",check.names = F,stringsAsFactors = F,header = T)
comgene <- intersect(rownames(Ginfo), rownames(tpm))
Ginfo <- Ginfo[comgene,]
tpm <- tpm[comgene,]
identical(rownames(tpm), rownames(Ginfo))
## [1] TRUE
tpm$Gene <- Ginfo[rownames(tpm),"gene"]
tpm <- as.data.frame(apply(tpm[,setdiff(colnames(tpm), "Gene")], 2, function(x) tapply(x, INDEX = factor(tpm$Gene), FUN=median, na.rm = TRUE))) # 对重复基因取表达谱中位数
colnames(tpm) <- substr(colnames(tpm),1,15)
expr <- tpm

# 提取突变和表达的共同样本
comsam <- intersect(colnames(expr), colnames(mut.binary))
mut.binary <- mut.binary[,comsam]
expr <- expr[,comsam]

# 加载原文用于聚类的差异表达基因
degs <- read.table("degs.txt",sep = "\t",row.names = NULL,header = T,check.names = F,stringsAsFactors = F)
comgene <- intersect(degs$Gene, rownames(expr))

# 从论文中获取iCluster亚型：Comprehensive and Integrative Genomic Characterization of Hepatocellular Carcinoma
icluster <- read.delim("icluster.txt",sep = "\t",row.names = 2,check.names = F,stringsAsFactors = F,header = T)
rownames(icluster) <- substr(rownames(icluster),1,15)

Nearest template prediction

加载不同HCC亚型签名并制作nearest template prediction所需的模版

boyault <- read.delim("Boyault-marker-all.txt",sep = "\t",row.names = NULL,check.names = F,stringsAsFactors = F,header = T)
tmp <- bitr(boyault$`Gene ID #1`,fromType = "ENTREZID",toType = "SYMBOL",OrgDb = org.Hs.eg.db)
## 'select()' returned 1:1 mapping between keys and columns
## Warning in bitr(boyault$`Gene ID #1`, fromType = "ENTREZID", toType =
## "SYMBOL", : 0.31% of input gene IDs are fail to map...
boyault <- boyault[,c(2,4)]
colnames(boyault) <- c("ENTREZID","class")
boyault <- merge(tmp,boyault,by = "ENTREZID", all.x = T)
colnames(boyault)[2] <- "probe"

chiang <- read.delim("Chiang-marker-all.txt",sep = "\t",row.names = NULL,check.names = F,stringsAsFactors = F,header = T)
chiang <- chiang[,c(1,4)]
colnames(chiang) <- c("probe","class")

hoshida <- read.delim("Hoshida-genelist.txt",sep = "\t",row.names = NULL,check.names = F,stringsAsFactors = F,header = T)
hoshida <- hoshida[,c(1,2)]
colnames(hoshida) <- c("probe","class")
hoshida$class <- paste0("C",hoshida$class)

lee <- read.delim("Lee-marker-all.txt",sep = "\t",row.names = NULL,check.names = F,stringsAsFactors = F,header = T)
lee <- lee[,c(1,3)]
colnames(lee) <- c("probe","class")

通过nearest template prediction获取其他亚型

tcga.boyault <- ntp(emat = t(scale(t(expr))),
                    templates = boyault,
                    doPlot = T,
                    seed = 19991018)

tcga.chiang <- ntp(emat = t(scale(t(expr))),
                   templates = chiang,
                   doPlot = T,
                   seed = 19991018)

tcga.lee <- ntp(emat = t(scale(t(expr))),
                   templates = lee,
                   doPlot = T,
                   seed = 19991018)

tcga.hoshida <- ntp(emat = t(scale(t(expr))),
                    templates = hoshida,
                    doPlot = T,
                    seed = 19991018)

开始画图

# 构建样本注释
annCol <- data.frame(CTNNB1 = ifelse(as.numeric(mut.binary["CTNNB1",]) == 0,"WT","MT"),
                     Lee = as.character(tcga.lee$prediction),
                     Hoshida = as.character(tcga.hoshida$prediction),
                     Boyault = as.character(tcga.boyault$prediction),
                     Chiang = as.character(tcga.chiang$prediction),
                     row.names = colnames(mut.binary),
                     stringsAsFactors = F)

annCol[intersect(rownames(annCol),rownames(icluster)),"iCluster"] <- icluster[intersect(rownames(annCol),rownames(icluster)), "iCluster clusters (k=3, Ronglai Shen)"]
annCol[is.na(annCol$iCluster),"iCluster"] <- "N/A"
annCol$Hoshida <- gsub("C","S",annCol$Hoshida)
annCol$iCluster <- gsub("iCluster:","iC",annCol$iCluster)
## 替换Chiang类型的标签
annCol[which(annCol$Chiang == "interferon class"), "Chiang"] <- "INTERFERON"
annCol[which(annCol$Chiang == "proliferation class"), "Chiang"] <- "PROLIFERATION"
annCol[which(annCol$Chiang == "CTNNB1 class"), "Chiang"] <- "CTNNB1"
annCol[which(annCol$Chiang == "unannotated class"), "Chiang"] <- "UNANNOTATED"
annCol[which(annCol$Chiang == "chromosome 7 polysomy class"), "Chiang"] <- "POLYSOMY7"

## 设置标签的因子顺序
annCol$CTNNB1 <- factor(annCol$CTNNB1, levels = c("MT","WT"))
annCol$Lee <- factor(annCol$Lee, levels = c("SURVIVAL_DN","SURVIVAL_UP"))
annCol$Hoshida <- factor(annCol$Hoshida, levels = c("S1","S2","S3"))
annCol$Boyault <- factor(annCol$Boyault, levels = c("G12","G3","G56"))
annCol$Chiang <- factor(annCol$Chiang, levels = c("PROLIFERATION","CTNNB1","INTERFERON","POLYSOMY7","UNANNOTATED"))
annCol$iCluster <- factor(annCol$iCluster, levels = c("iC1","iC2","iC3","N/A"))

# 构建样本注释所对应的颜色
annColors <- list()
annColors[["CTNNB1"]] <- c("WT" = "white","MT" = "black")
annColors[["Lee"]] <- c("SURVIVAL_DN" = "#E80035","SURVIVAL_UP" = "#DBDEDD")
annColors[["Hoshida"]] <- c("S1" = "#E8536B","S2" = "#F6B879","S3" = "#DFEBAF")
annColors[["Boyault"]] <- c("G12" = "#E8536B","G3" = "#F6B879","G56" = "#DFEBAF")
annColors[["Chiang"]] <- c("PROLIFERATION" = "#F2A1A2","CTNNB1" = "#68BE8B","INTERFERON" = "#2CA8E1","POLYSOMY7" = "#A2D9F1","UNANNOTATED" = "#BAC8E5")
annColors[["iCluster"]] <- c("iC1" = "#A8BFDE","iC2" = "#DA8A88","iC3" = "#BCA3D6","N/A" = "white")
annColors[["Group"]] <- c("A" = "#FE0000","B" = "black")
annColors[["MS"]] <- c("MS1" = "#FE0000", "MS2" = "#00AF50","MS3" = "#0071C0")

# 对表达谱进行无监督聚类
indata <- expr[comgene,]
indata <- indata[rowSums(indata) > 0,]
hcs <- hclust(distanceMatrix(as.matrix(indata), "euclidean"), "ward.D")
hcg <- hclust(distanceMatrix(as.matrix(t(indata)), "euclidean"), "ward.D")

group <- cutree(hcs, k = 3)
table(group, annCol$CTNNB1) # 查看一下样本量以及CTNNB1突变情况来确定以下的类型命名
##      
## group MT WT
##     1 47  0
##     2 32 61
##     3 15 28
annCol$MS <- ifelse(group == 1, "MS2", ifelse(group == 3, "MS1","MS3")) # 由于1组有更多CTNNB1，所以是MS2，由于3组样本更少，所以是MS1
annCol$Group <- ifelse(group == 3, "A", "B") # 由于3组是MS1所以为Group A

plotdata <- standarize.fun(indata, halfwidth = 2) # 数据标准化用于绘图
hm <- pheatmap(plotdata,
         color = bluered(64),
         border_color = NA,
         cluster_rows = hcg,
         cluster_cols = hcs,
         treeheight_row = 30,
         treeheight_col = 30,
         #cutree_cols = 3,
         show_rownames = F,
         show_colnames = F,
         annotation_col = annCol[,c("iCluster","Chiang","Boyault","Hoshida","Lee","CTNNB1","MS","Group")],
         annotation_colors = annColors,
         cellwidth = 0.6,
         cellheight = 0.2)

pdf("heatmap.pdf", width = 8,height = 6)
draw(hm, heatmap_legend_side = "left", annotation_legend_side = "left")
invisible(dev.off())

计算group和其他类型的独立性

outTabGroup <- outTabMS <- NULL
for (i in c("iCluster","Chiang","Boyault","Hoshida","Lee","CTNNB1")) {
  tmp <- annCol[,c(i,"Group")]
  set.seed(123)
  p <- fisher.test(table(tmp[,1],tmp[,2]), 
                   simulate.p.value = T)$p.value # 由于可能存在比较复杂的列联表，所以这里直接使用模拟的p值
  outTabGroup <- rbind.data.frame(outTabGroup,
                                  data.frame(VarA = "Group",
                                             VarB = i,
                                             p = p),
                                  stringsAsFactors = F)
  
  tmp <- annCol[,c(i,"MS")]
  set.seed(123)
  p <- fisher.test(table(tmp[,1],tmp[,2]),simulate.p.value = T)$p.value
  outTabMS <- rbind.data.frame(outTabMS,
                               data.frame(VarA = "MS",
                                          VarB = i,
                                          p = p),
                               stringsAsFactors = F)
}

# 把p value保存到文件
write.table(outTabGroup, file = "independent test between group and other classification.txt",sep = "\t",row.names = F,col.names = T,quote = F)
write.table(outTabMS, file = "independent test between ms and other classification.txt",sep = "\t",row.names = F,col.names = T,quote = F)

# 保存镜像
#save.image(file = "LIHC.RData")

后期加工

输出的PDF文件是矢量图，可以用Illustrator等矢量图编辑工具打开，添加p value。如果想用代码添加p value，可参考FigureYa165heatmapPvalue或FigureYa280TMEofSTS

Session Info

sessionInfo()
## R version 4.1.1 (2021-08-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] gplots_3.1.3          ComplexHeatmap_2.8.0  ClassDiscovery_3.4.0 
##  [4] oompaBase_3.2.9       cluster_2.1.3         org.Hs.eg.db_3.13.0  
##  [7] AnnotationDbi_1.54.1  IRanges_2.26.0        S4Vectors_0.30.2     
## [10] Biobase_2.54.0        BiocGenerics_0.40.0   clusterProfiler_4.0.5
## [13] forcats_0.5.1         stringr_1.4.0         dplyr_1.0.9          
## [16] purrr_0.3.4           readr_2.1.2           tidyr_1.2.0          
## [19] tibble_3.1.8          ggplot2_3.3.6         tidyverse_1.3.2      
## [22] CMScaller_2.0.1      
## 
## loaded via a namespace (and not attached):
##   [1] circlize_0.4.15        readxl_1.4.0           shadowtext_0.1.2      
##   [4] backports_1.4.1        fastmatch_1.1-3        plyr_1.8.7            
##   [7] igraph_1.3.4           lazyeval_0.2.2         splines_4.1.1         
##  [10] BiocParallel_1.26.2    GenomeInfoDb_1.28.4    digest_0.6.29         
##  [13] foreach_1.5.2          yulab.utils_0.0.5      htmltools_0.5.3       
##  [16] GOSemSim_2.18.1        viridis_0.6.2          GO.db_3.13.0          
##  [19] fansi_1.0.3            magrittr_2.0.3         memoise_2.0.1         
##  [22] doParallel_1.0.17      googlesheets4_1.0.0    tzdb_0.3.0            
##  [25] Biostrings_2.60.2      graphlayouts_0.8.0     modelr_0.1.8          
##  [28] matrixStats_0.62.0     vroom_1.5.7            enrichplot_1.12.3     
##  [31] colorspace_2.0-3       blob_1.2.3             rvest_1.0.2           
##  [34] ggrepel_0.9.1          haven_2.5.0            xfun_0.31             
##  [37] crayon_1.5.1           RCurl_1.98-1.8         jsonlite_1.8.0        
##  [40] scatterpie_0.1.7       iterators_1.0.14       ape_5.6-2             
##  [43] glue_1.6.2             polyclip_1.10-0        gtable_0.3.0          
##  [46] gargle_1.2.0           zlibbioc_1.38.0        XVector_0.32.0        
##  [49] GetoptLong_1.0.5       shape_1.4.6            scales_1.2.0          
##  [52] DOSE_3.18.3            DBI_1.1.3              Rcpp_1.0.9            
##  [55] viridisLite_0.4.0      clue_0.3-61            gridGraphics_0.5-1    
##  [58] tidytree_0.3.9         mclust_5.4.10          bit_4.0.4             
##  [61] httr_1.4.3             fgsea_1.18.0           RColorBrewer_1.1-3    
##  [64] ellipsis_0.3.2         pkgconfig_2.0.3        farver_2.1.1          
##  [67] sass_0.4.2             dbplyr_2.2.1           utf8_1.2.2            
##  [70] ggplotify_0.1.0        tidyselect_1.1.2       rlang_1.0.4           
##  [73] reshape2_1.4.4         munsell_0.5.0          cellranger_1.1.0      
##  [76] tools_4.1.1            cachem_1.0.6           oompaData_3.1.2       
##  [79] downloader_0.4         cli_3.3.0              generics_0.1.3        
##  [82] RSQLite_2.2.15         broom_1.0.0            evaluate_0.16         
##  [85] fastmap_1.1.0          yaml_2.3.5             ggtree_3.0.4          
##  [88] knitr_1.39             bit64_4.0.5            fs_1.5.2              
##  [91] tidygraph_1.2.1        caTools_1.18.2         KEGGREST_1.32.0       
##  [94] ggraph_2.0.6           nlme_3.1-159           aplot_0.1.6           
##  [97] DO.db_2.9              xml2_1.3.3             compiler_4.1.1        
## [100] rstudioapi_0.13        png_0.1-7              reprex_2.0.1          
## [103] treeio_1.16.2          tweenr_1.0.2           bslib_0.4.0           
## [106] stringi_1.7.8          highr_0.9              lattice_0.20-45       
## [109] Matrix_1.4-1           vctrs_0.4.1            pillar_1.8.0          
## [112] lifecycle_1.0.1        GlobalOptions_0.1.2    jquerylib_0.1.4       
## [115] data.table_1.14.2      cowplot_1.1.1          bitops_1.0-7          
## [118] patchwork_1.1.1        qvalue_2.24.0          R6_2.5.1              
## [121] KernSmooth_2.23-20     gridExtra_2.3          codetools_0.2-18      
## [124] gtools_3.9.3           MASS_7.3-58.1          assertthat_0.2.1      
## [127] rjson_0.2.21           withr_2.5.0            GenomeInfoDbData_1.2.6
## [130] parallel_4.1.1         hms_1.1.1              ggfun_0.0.6           
## [133] rmarkdown_2.14         googledrive_2.0.0      Cairo_1.6-0           
## [136] ggforce_0.3.3          lubridate_1.8.0